Mapeamento digital do solo para agricultura de precisão
FEAGRI-UNICAMP
Alessandro Samuel-Rosa
25 e 26 de janeiro de 2019
Este estudo tem como objetivo realizar o mapeamento de uma série de propriedades do solo em uma fazenda de produção de bovinos de corte no município de Caiuá, no oeste do estado de São Paulo. Dentre as propriedades do solo estão: conteúdo de fósforo e bases trocáveis, capacidade de troca de cátions e saturação por bases, conteúdo de micronutrientes, pH e acidez potencial, conteúdo de carbono orgânico, e distribuição do tamanho de partículas.
1 Material e Métodos
1.1 Área de estudo
O primeiro passo consiste em fazer um levantamento das informações disponíveis sobre a área de estudo. Essas informações incluem dados espacialmente exaustivos, como imagens de satélite e índices de vegetação, modelos digitais de elevação e atributos de terreno, mapas pedológicos e geológicos, entre outros. Também incluem informações históricas sobre o uso da terra e as práticas de manejo e conservação do solo.
Vamos iniciar visualizando a área de estudo sobre uma imagem de satélite da coleção pública da ESRI. Para isso, carregamos o arquivo vetorial contendo o limite da área de estudo, "../data/vector/farm.shp", usando a função sf::read_sf. Em seguida, criamos uma figura responsiva usando a função mapview::mapview.
farm <- read_sf("../data/vector/farm.shp")
mapviewFarm <-
function () {
mapview(farm, color = "red", alpha.regions = 0.01, legend = FALSE)
}
mapviewFarm()1.1.1 Relevo
A imagem acima mostra que a área de estudo é cortada por inúmeros terraços, estruturas construídas para reduzir a energia cinética da água da chuva escoando na superfície do solo e, assim, reduzir a erosão do solo. Também observamos algumas escavações para coleta da enxurrada em pontos de convergência do terreno. Como os fluxos hídricos na paisagem influenciam a formação do solo, espera-se que o relevo da área de estudo exerça algum controle sobre a distribuição espacial das propriedades do solo.
Vamos analisar mais de perto o relevo da área de estudo utilizando o modelo digital de elevação (MDE). Para isso, primeiro calculamos o sombreamento do relevo usando sf::gdal_utils. Em seguida, carregamos os dois arquivos, "../data/raster/elevation.tif" e "../data/raster/hillshade.tif", usando raster::stack, para então gerar uma imagem com as duas camadas sobrepostas, dando assim a impressão visual do formato tridimensional do terreno. Para auxiliar, adicionamos curvas de nível geradas com a função raster::rasterToContour.
gdal_utils(
util = "demprocessing", source = "../data/raster/elevation.tif", destination = "../data/raster/hillshade.tif",
processing = "hillshade")
dem <- stack("../data/raster/elevation.tif", "../data/raster/hillshade.tif")
contour_lines <- dem[["elevation"]] %>% rasterToContour() %>% st_as_sf()
plot(dem[["hillshade"]], col = grey(0:100/100), legend = FALSE)
plot(dem[["elevation"]], col = topo.colors(24, alpha = 0.50), add = TRUE)
plot(contour_lines, add = TRUE, col = "black")
Exercício 1 A imagem acima mostra a configuração do terreno na área de estudo. Você espera que o terreno esteja relacionado à distribuição espacial das propriedades do solo que se pretende mapear? De que maneira se daria essa relação? Quais atributos de terreno podem ser úteis para descrever o efeito do terreno sobre a distribuição espacial das propriedades do solo? Faça uma lista com, no mínimo, cinco atributos de terreno.
- …
- …
- …
- …
- …
Atributos de terreno podem ser computados usando o SAGA GIS. Por exemplo, RSAGA::rsaga.slope.asp.curv calcula alguns atributos locais do terreno, como a declividade e a curvatura do terreno. Alguns atributos regionais do terreno, como a área de captação e o índice de umidade topográfica, podem ser calculados usando RSAGA::rsaga.wetness.index. Note que o formato de arquivo matricial do SAGA GIS é SGRD – para escrita – e SDAT – para leitura.
rsaga.slope.asp.curv(
in.dem = "../data/raster/elevation.tif",
out.slope = "../data/raster/slope.sgrd",
out.cgene = "../data/raster/general_curvature.sgrd",
out.cplan = "../data/raster/plan_curvature.sgrd",
out.cprof = "../data/raster/profile_curvature.sgrd",
unit.slope = "percent")## Search for SAGA command line program and modules...
## Done
##
##
## SAGA Version: 2.3.1
##
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_morphometry
## library : Morphometry
## tool : Slope, Aspect, Curvature
## author : O.Conrad (c) 2001
## processors : 4 [4]
##
## loading: elevation
##
##
## Parameters
##
##
## Grid system: 29.286570; 81x 38y; 404060.656585x 7606472.659017y
## Elevation: elevation
## Slope: Slope
## Aspect: Aspect
## General Curvature: General Curvature
## Profile Curvature: Profile Curvature
## Plan Curvature: Plan Curvature
## Tangential Curvature: <not set>
## Longitudinal Curvature: <not set>
## Cross-Sectional Curvature: <not set>
## Minimal Curvature: <not set>
## Maximal Curvature: <not set>
## Total Curvature: <not set>
## Flow Line Curvature: <not set>
## Method: 9 parameter 2nd order polynom (Zevenbergen & Thorne 1987)
## Slope Units: percent
## Aspect Units: radians
##
## Save grid: ../data/raster/slope.sgrd...
## Save grid: /tmp/RtmpWszgZ3/file128115a11c3d...
## Save grid: ../data/raster/general_curvature.sgrd...
## Save grid: ../data/raster/profile_curvature.sgrd...
## Save grid: ../data/raster/plan_curvature.sgrd...## Warning in rsaga.slope.asp.curv(in.dem = "../data/raster/elevation.tif", : Plan and profile curvature calculations have changed with SAGA 2.1.1+
## See help(rsaga.slope.asp.curv) for more informationrsaga.wetness.index(
in.dem = "../data/raster/elevation.tif",
out.wetness.index = "../data/raster/twi.sgrd",
out.carea = "../data/raster/catchment_area.sgrd",
out.cslope = "../data/raster/catchment_slope.sgrd",
area.type = "absolute")## Search for SAGA command line program and modules...
## Done
##
##
## SAGA Version: 2.3.1
##
## library path: /usr/lib/x86_64-linux-gnu/saga/
## library name: libta_hydrology
## library : Hydrology
## tool : SAGA Wetness Index
## author : (c) 2001 by J.Boehner, O.Conrad
## processors : 4 [4]
##
## loading: elevation
##
##
## Parameters
##
##
## Grid system: 29.286570; 81x 38y; 404060.656585x 7606472.659017y
## Elevation: elevation
## Weights: <not set>
## Catchment area: Catchment area
## Catchment slope: Catchment slope
## Modified Catchment Area: Modified Catchment Area
## Topographic Wetness Index: Topographic Wetness Index
## Suction: 10.000000
## Type of Area: absolute catchment area
## Type of Slope: catchment slope
## Minimum Slope: 0.000000
## Offset Slope: 0.100000
## Slope Weighting: 1.000000
##
## Create index: elevation
## catchment area and slope...
## pass 1 (1560 > 0)
## pass 2 (1203 > 0)
## pass 3 (764 > 0)
## pass 4 (502 > 0)
## pass 5 (381 > 0)
## pass 6 (230 > 0)
## pass 7 (149 > 0)
## pass 8 (94 > 0)
## pass 9 (57 > 0)
## pass 10 (40 > 0)
## pass 11 (30 > 0)
## pass 12 (19 > 0)
## pass 13 (14 > 0)
## pass 14 (1 > 0)
## pass 15 (0 > 0)
## post-processing...
## topographic wetness index...
## Save grid: ../data/raster/catchment_area.sgrd...
## Save grid: ../data/raster/catchment_slope.sgrd...
## Save grid: /tmp/RtmpWszgZ3/file128185cb4bf...
## Save grid: ../data/raster/twi.sgrd...Em seguida carregamos os atributos de terreno computados usando o SAGA GIS para dentro do objeto dem – criado acima – usando raster::stack. Calculamos o logaritmo da área de captação para transformar os dados para uma escala mais estreita e interpretável. Ainda, transformamos os valores de declividade da área de captação de radianos para porcentagem.
dem <-
stack(
c("../data/raster/slope.sdat",
"../data/raster/general_curvature.sdat",
"../data/raster/plan_curvature.sdat",
"../data/raster/profile_curvature.sdat",
"../data/raster/twi.sdat",
"../data/raster/catchment_area.sdat",
"../data/raster/catchment_slope.sdat")) %>%
stack(dem)
dem[["catchment_area"]] <- log1p(dem[["catchment_area"]])
dem[["catchment_slope"]] <- tan(dem[["catchment_slope"]]) * 100Depois de processados, podemos visualizar os atributos do terreno usando o sombreamento e as curvas de nível para auxiliar na sua interpretação. Como são vários os atributos do terreno, precisamos iterar sobre a função plot tantas vezes quantos forem os atributos do terreno usando a função for.
nr <- nlayers(dem) %>% sqrt() %>% ceiling()
nc <- nlayers(dem) / nr
par(mfrow = c(nr, nc))
for (i in nlayers(dem):1) {
plot(dem[["hillshade"]], col = grey(0:100/100), legend = FALSE, main = names(dem)[i], axes = FALSE)
plot(dem[[i]], col = terrain.colors(24, alpha = 0.50), add = TRUE)
plot(contour_lines, add = TRUE, col = "black", reset = FALSE)
}
1.1.2 Uso da terra e manejo do solo
Vamos retornar à imagem de satélite da área de estudo para identificar feições relacionadas ao uso da terra e às práticas de manejo e conservação do solo que podem exercer algum controle sobre a distribuição espacial das propriedades do solo. Adicionamos as curvas de nível para auxiliar na representação da forma do terreno e, assim, identificar possíveis correlações entre o relevo e o uso da terra e as práticas de manejo e conservação do solo.
mapview(contour_lines) +
mapviewFarm()Exercício 2 A imagem de satélite nos permite identificar uma série de feições espaciais relacionadas ao uso da terra e às práticas de manejo e conservação do solo na área de estudo. Quais dessas feições você espera que estejam relacionadas aos valores – a serem observados – das propriedades do solo que se pretende mapear e, assim, estar determinando a sua distribuição espacial? Faça uma lista com, no mínimo, 10 itens.
- …
- …
- …
- …
- …
- …
- …
- …
- …
- …
A primeira feição espacial relacionada ao uso da terra identificada na área de estudo é a sua subdivisão em piquetes de pastoreio, indicando que os animais são conduzidos num sistema de pastoreio rotativo. O formato retangular dos piquetes não indica nenhuma relação entre a sua construção e as feições do terreno. Contudo, é possível que diferentes práticas de manejo e conservação do solo sejam aplicadas em cada piquete, influenciando assim a distribuição espacial das propriedades do solo.
Vamos carregar o arquivo vetorial ../data/vector/fields.shp contendo os limites dos 26 piquetes de pastoreio e as quatro áreas de circulação interna dos animais. Em seguida, visualizamos os piquetes sobre a imagem de satélite usando mapview::mapview. Adicionamos o limite da área de estudo para verificar a consistência espacial dos dados dos piquetes de pastoreio.
fields <- read_sf("../data/vector/fields.shp")
mapview(fields) +
mapviewFarm()Para utilizar as informações sobre os piquetes nas etapas subsequentes, é preciso transformar os dados vetoriais em dados matriciais. Vamos usar a função raster::rasterize para isso, usando o objeto dem como referência.
fields %<>%
rasterize(y = dem[["slope"]], field = "id")
names(fields) <- "fields"
mapview(fields) +
mapviewFarm()Uma segunda informação que pode importante para explicar a distribuição espacial das propriedades do solo é a localização dos bebedouros. Isso porque os bebedouros, e seu entorno, são locais onde os animais costumam se aglomerar periodicamente. Espera-se que pisoteio e a defecação e micção periódicos nesses locais altere, consideravelmente, as propriedades do solo, criando um ambiente diferenciado em relação ao restante da área de estudo.
Vamos criar um vetor de pontos para indicar a localização dos tanques de água usando a função mapedit::drawFeatures. Note que mapedit::drawFeatures retorna um objeto com coordenadas geográficas. Para transformar as coordenadas para o sistema de referência de coodernadas EPSG:32722, usamos a função sf::st_transform. Ainda, criamos uma coluna id para identificar cada um dos tanques (dplyr::mutate) e descartamos as colunas criadas por mapedit::drawFeatures (dplyr::select). Salvamos o arquivo vetorial em disco – "../data/vector/water_tank.shp" para facilitar o seu reúso.
if (!file.exists("../data/vector/water_tank.shp")) {
water_tank <-
mapview(farm, alpha.regions = 0.01, color = "red") %>%
drawFeatures() %>%
st_transform(crs = 32722) %>%
dplyr::mutate(id = 1:nrow(.)) %>%
dplyr::select(id)
write_sf(water_tank, "../data/vector/water_tank.shp")
mapview(water_tank) +
mapviewFarm()
} else {
water_tank <- read_sf("../data/vector/water_tank.shp")
mapview(water_tank) +
mapviewFarm()
}Assim como fora feito com os polígonos que representam os piquetes de pastoreio, os dados sobre a localização dos bebedouros precisam ser transformados para o formato matricial. Como se trata de dados pontuais, que não cobrem toda a área de estudo, uma solução é computar uma mapa de distâncias até os tanques de água. Isso é razoável porque, em princípio, espera-se que o efeito do pisoteio, defecação e micção seja bastante localizado, diminuindo a medida que nos afastamos dos bebedouros. Para calcular a distância até os bebedouros usamos a função raster::distanceFromPoints – note que, primeiro, é preciso converter a classe do objeto water_tank de sf para SpatialPointsDataFrame.
water_tank_dist <-
water_tank %>%
as_Spatial() %>%
distanceFromPoints(dem[["slope"]], .)
names(water_tank_dist) <- "water_tank_dist"
mapview(water_tank_dist) +
mapview(water_tank, legend = FALSE) +
mapviewFarm()1.1.3 Imagens de satélite e índices de vegetação
A análise da imagens de satélite da área de estudo mostra que a pastagem exibe alguma variação espacial na sua coloração. Essa variação pode ser devida às propriedades do solo. Assim, as bandas de um sensor orbital e os índices de vegetação podem ser úteis para explicar a variação espacial das propriedades do solo. Está disponível o índice de vegetação por diferença normalizada (NDVI) da área de estudo em três datas: 22/04/2013, 25/01/2016 e 20/04/2018. Vamos carregar os dados usando raster::stack.
ndvi <-
c("2013_04_22", "2016_01_25", "2018_04_20") %>%
paste("../data/raster/ndvi_", ., ".tif", sep = "") %>%
stack() %T>%
plot()
A distribuição espacial do NDVI parece bastante correlacionada com a disposição dos piquetes de pastoreio. Contudo, seria interessante conhecer a variação temporal do NDVI. Para isso, podemos calcular o seu desvio padrão em cada célula usando a função raster::calc. Quanto maior o desvio padrão, maior a variação temporal do NDVI. A magnitude da variação temporal do NDVI pode ajudar a identificar, por exemplo, áreas onde o solo impõe maiores restrições ao crescimento vegetativo ao longo dos anos como menor fertilidade, menor profundidade, menor capacidade de retenção de água e elevada pedregosidade. Em locais como esses, a vegetação costuma ser mais sensível às condições meteorológicas extremas, apresentando mais variação no NDVI. (Note que os dados estão disponíveis para apenas três anos, o que pode não ser suficiente para representar a variação temporal do NDVI de maneira precisa.)
ndvi$ndvi_variation <-
ndvi %>%
calc(sd) %T>%
plot()
1.1.4 Propriedades do solo
Informações indiretas sobre as propriedades do solo também podem auxiliar no entendimento da variação espacial das propriedades do solo que se pretende mapear. Uma delas é a condutividade elétrica aparente do solo. Vamos carregar os dados de condutividade elétrica aparente – "../data/vector/ec.shp" – disponíveis para a área de estudo. Note que os dados estão no formato vetorial (pontos) e precisam ser transformados para o formato matricial via interpolação. Contudo, como a densidade de pontos é muito elevada, faremos uma filtragem dos dados antes de proceder com a interpolação propriamente dita – com raster::rasterize. Para isso, usaremos a mediana dos valores observados dentro de cada célula de 30 m x 30 m, ignorando assim valores discrepantes. Em seguida, procedemos com a interpolação determinística usando as funções fields::Tps (thin plate spline) e raster::interpolate.
if (!file.exists("../data/raster/ec.grd")) {
ec <-
read_sf("../data/vector/ec.shp") %>%
dplyr::select(CE37, CE75) %>%
rasterize(dem[["slope"]], fun = median) %>%
as("SpatialPointsDataFrame")
tps_ec37 <- Tps(ec@coords, ec$CE37)
tps_ec75 <- Tps(ec@coords, ec$CE75)
ec %<>% rasterize(dem[["slope"]])
ec <- ec[[c("CE37", "CE75")]]
ec[["CE37"]] <- interpolate(ec[["CE37"]], tps_ec37)
ec[["CE75"]] <- interpolate(ec[["CE75"]], tps_ec75)
writeRaster(ec, file = "../data/raster/ec.grd", format = "raster")
} else {
ec <- stack("../data/raster/ec.grd")
}
plot(ec)
1.1.5 Covariáveis
Para finalizar o processamento das covariáveis, vamos juntá-las em um único objeto covar_data. Os três objetos originais, dem, ec e ndvi, podem ser removidos do ambiete de trabalho e assim liberar espaço da memória. Finalmente, adicionamos as coordenadas espaciais, x e y, às covariáveis a fim de que também sejam usadas para calibração dos modelos e predição espacial das propriedades do solo.
covar_data <-
stack(dem, ec, ndvi, fields, water_tank_dist) %>%
dropLayer("hillshade")
covar_data$x <- xFromCell(covar_data, 1:ncell(covar_data))
covar_data$y <- yFromCell(covar_data, 1:ncell(covar_data))
rm(dem, ec, ndvi, fields, water_tank_dist)A correlação linear entre as covariáveis pode ser observada usando a função pedometrics::plotCor. Note que, em geral, a correlação linear entre as covariáveis é pequena, o que indica que há pouca redundância nos dados. As maiores correlações observada são entre a curvatura geral e a curvatura de perfil (r = 0.90) e entre a elevação e a coordenada y (r = -0.96).
covar_data %>%
values() %>%
cor(use = "complete") %>%
round(2) %>%
plotCor(cex = 0.75)## Too long column names found. Replacing with x1 - slope; x2 - general_curvature; x3 - plan_curvature; x4 - profile_curvature; x5 - twi; x6 - catchment_area; x7 - catchment_slope; x8 - elevation; x9 - CE37; x10 - CE75; x11 - ndvi_2013_04_22; x12 - ndvi_2016_01_25; x13 - ndvi_2018_04_20; x14 - ndvi_variation; x15 - fields; x16 - water_tank_dist; x17 - x; x18 - y.
1.2 Amostragem do solo
Uma vez construída a base de dados espaciais da área de estudo, podemos proceder com a escolha dos locais de observação do solo. Isso pode ser feito usando algoritmos que escolhem os locais de observação do solo de maneira a atender a algum critério objetivo. O pacote spsann possui uma série de funções com essa finalidade. Uma delas, spsann::optimDIST, otimiza uma configuração amostral de maneira a melhor representar a distribuição marginal empírica de cada uma das covariáveis. Se entendemos que as covariáveis estão intimamente relacionadas à distribuição espacial das propriedades do solo, então, ao representarmos sua distribuição marginal empírica, aumentamos nossas chances de obter dados mais representativos e, assim, produzir mapas mais acurados das propriedades do solo.
Os dados do solo utilizados neste estudo já foram coletados. Portanto, deixaremos de lado o processo de otimização e passaremos imediatamente ao processo de calibração dos modelos, validação e predição espacial.
Os dados do solo estão numa planilha no Google Drive (https://goo.gl/qcM2uq), organizados conforme os padrões do Repositório Brasileiro Livre para Dados Abertos do Solo (http://coral.ufsm.br/febr/). Para descarregá-los, usaremos as funções googlesheets::gs_key e googlesheets::gs_read_csv. Note que é preciso usar a função readr::locale para especificar que utilizamos a vírgula como separador decimal.
São duas as planilhas com dados, observacao e camada, que fundimos usando merge e, em seguida, transformamos para a classe sf usando sf::st_as_sf.
locale <- locale(decimal_mark = ",")
camada <-
"1VqP_W9rS4DJN3EwrGnkT_gFduMTwaupFgct-9OqK9WI" %>%
gs_key() %>%
gs_read_csv(ws = "camada", locale = locale, comment = "#metadado>", verbose = FALSE)## Warning: Duplicated column names deduplicated: 'indefinido' =>
## 'indefinido_1' [29], 'indefinido' => 'indefinido_2' [30], 'indefinido' =>
## 'indefinido_3' [31], 'indefinido' => 'indefinido_4' [32], 'indefinido' =>
## 'indefinido_5' [33], 'indefinido' => 'indefinido_6' [34], 'indefinido' =>
## 'indefinido_7' [35], 'indefinido' => 'indefinido_8' [36], 'indefinido' =>
## 'indefinido_9' [37]observacao <-
"1VqP_W9rS4DJN3EwrGnkT_gFduMTwaupFgct-9OqK9WI" %>%
gs_key() %>%
gs_read_csv(ws = "observacao", locale = locale, comment = "#metadado>", verbose = FALSE)
soil_data <-
merge(observacao, camada, by = "observacao_id") %>%
st_as_sf(coords = c("coord_x", "coord_y"), crs = 32722) %T>%
print()## Simple feature collection with 627 features and 50 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 404071 ymin: 7606493 xmax: 406374 ymax: 7607536
## epsg (SRID): 32722
## proj4string: +proj=utm +zone=22 +south +datum=WGS84 +units=m +no_defs
## First 10 features:
## observacao_id observacao_data coord_sistema coord_precisao coord_fonte
## 1 8 xx-xx-2018 EPSG:32722 NA GPS
## 2 8 xx-xx-2018 EPSG:32722 NA GPS
## 3 35 xx-xx-2018 EPSG:32722 NA GPS
## 4 35 xx-xx-2018 EPSG:32722 NA GPS
## 5 48 xx-xx-2018 EPSG:32722 NA GPS
## 6 48 xx-xx-2018 EPSG:32722 NA GPS
## 7 54 xx-xx-2018 EPSG:32722 NA GPS
## 8 54 xx-xx-2018 EPSG:32722 NA GPS
## 9 57 xx-xx-2018 EPSG:32722 NA GPS
## 10 65 xx-xx-2018 EPSG:32722 NA GPS
## pais_id estado_id municipio_id amostra_tipo amostra_quanti amostra_area
## 1 BR SP Caiuá COMPOSTA 4 78.54
## 2 BR SP Caiuá COMPOSTA 4 78.54
## 3 BR SP Caiuá COMPOSTA 4 78.54
## 4 BR SP Caiuá COMPOSTA 4 78.54
## 5 BR SP Caiuá COMPOSTA 4 78.54
## 6 BR SP Caiuá COMPOSTA 4 78.54
## 7 BR SP Caiuá COMPOSTA 4 78.54
## 8 BR SP Caiuá COMPOSTA 4 78.54
## 9 BR SP Caiuá COMPOSTA 4 78.54
## 10 BR SP Caiuá COMPOSTA 4 78.54
## camada_id amostra_id camada_nome profund_sup profund_inf
## 1 1 277648/2018 - 0 20
## 2 2 277649/2018 - 20 40
## 3 1 277650/2018 - 0 20
## 4 2 277651/2018 - 20 40
## 5 1 277652/2018 - 0 20
## 6 2 277653/2018 - 20 40
## 7 1 277654/2018 - 0 20
## 8 2 277655/2018 - 20 40
## 9 1 277656/2018 - 0 20
## 10 1 277657/2018 - 0 20
## fosforo_resina matorg_cromo carbono_xxx_xxx_xxx ph_cacl2 ph_smp
## 1 4 14 8 4.9 6.95
## 2 2 10 6 4.6 6.95
## 3 7 16 9 5.1 7.04
## 4 3 11 6 5.0 6.99
## 5 5 15 9 5.0 6.90
## 6 2 12 7 5.1 7.17
## 7 4 16 9 5.0 7.03
## 8 2 11 6 4.9 7.11
## 9 4 14 8 5.0 6.98
## 10 6 12 7 5.2 7.21
## potassio_resina calcio_resina magnesio_resina sodio_resina
## 1 1.1 13 8 0.2
## 2 0.6 15 6 0.3
## 3 2.7 22 9 0.3
## 4 0.7 25 6 0.2
## 5 0.8 23 8 0.1
## 6 0.6 26 6 0.2
## 7 1.0 23 8 0.2
## 8 0.9 18 7 0.1
## 9 1.2 14 9 0.2
## 10 2.8 15 9 0.3
## acidez_xxx_xxx aluminio_xxx_xxx indefinido ctc_soma_calc
## 1 16 0 16 38.3
## 2 16 0 16 37.9
## 3 14 0 14 48.0
## 4 15 0 15 46.9
## 5 16 0 16 47.9
## 6 12 0 12 44.8
## 7 14 0 14 46.2
## 8 13 0 13 39.0
## 9 15 0 15 39.4
## 10 12 0 12 39.1
## bases_soma_calc bases_saturacao_calc aluminio_saturacao_calc
## 1 22.3 58 0
## 2 21.9 58 0
## 3 34.0 71 0
## 4 31.9 68 0
## 5 31.9 67 0
## 6 32.8 73 0
## 7 32.2 70 0
## 8 26.0 67 0
## 9 24.4 62 0
## 10 27.1 69 0
## enxofre_fosfca boro_h2o cobre_dtpa_eaa ferro_dtpa_eaa manganes_dtpa_eaa
## 1 12 0.50 1.1 7 16.9
## 2 18 0.85 1.5 11 12.4
## 3 13 0.79 1.7 12 31.9
## 4 16 0.56 2.4 10 30.3
## 5 19 0.65 1.1 23 23.0
## 6 17 0.68 1.4 10 18.0
## 7 12 0.79 1.1 18 23.1
## 8 16 0.67 1.1 11 24.3
## 9 15 0.59 1.6 17 24.4
## 10 17 0.64 1.6 9 20.1
## zinco_dtpa_eaa indefinido_1 indefinido_2 indefinido_3 indefinido_4
## 1 0.8 2.9 33.9 20.9 0.5
## 2 0.4 1.6 39.6 15.8 0.8
## 3 2.0 5.6 45.8 18.8 0.6
## 4 2.4 1.5 53.3 12.8 0.4
## 5 1.2 1.7 48.0 16.7 0.2
## 6 0.7 1.3 58.0 13.4 0.4
## 7 2.1 2.2 49.8 17.3 0.4
## 8 1.4 2.3 46.2 17.9 0.3
## 9 2.9 3.0 35.5 22.8 0.5
## 10 1.8 7.2 38.4 23.0 0.8
## indefinido_5 indefinido_6 indefinido_7 indefinido_8 indefinido_9
## 1 0 41.8 11.8 1.6 7.3
## 2 0 42.2 25.0 2.5 10.0
## 3 0 29.2 8.1 2.4 3.3
## 4 0 32.0 35.7 4.2 8.6
## 5 0 33.4 28.8 2.9 10.0
## 6 0 26.8 43.3 4.3 10.0
## 7 0 30.3 23.0 2.9 8.0
## 8 0 33.3 20.0 2.6 7.8
## 9 0 38.1 11.7 1.6 7.5
## 10 0 30.7 5.4 1.7 3.2
## argila_xxx_xxx_xxx silte_xxx_xxx_xxx areia_xxx_xxx_xxx
## 1 150 121 729
## 2 199 94 707
## 3 136 116 748
## 4 164 77 759
## 5 143 63 794
## 6 179 96 725
## 7 106 130 764
## 8 148 85 767
## 9 135 73 792
## 10 111 76 813
## geometry
## 1 POINT (404358.2 7607536)
## 2 POINT (404358.2 7607536)
## 3 POINT (404568.9 7607508)
## 4 POINT (404568.9 7607508)
## 5 POINT (404310 7607486)
## 6 POINT (404310 7607486)
## 7 POINT (404433.9 7607499)
## 8 POINT (404433.9 7607499)
## 9 POINT (404495.3 7607501)
## 10 POINT (404654.7 7607490)Vamos visualizar a distribuição espacial das observações do solo na área de estudo. São 399 amostras da camada de 0-20 cm e 228 amostras da camada de 20-40 cm.
n <-
soil_data %>%
mutate(camada_id = as.factor(camada_id)) %>%
summarise(summary(camada_id)) %>%
unlist()
mapviewFarm() +
mapview(soil_data[soil_data$camada_id == 1, "observacao_id"], layer.name = glue("0-20 cm<br>n = {n[1]}")) +
mapview(soil_data[soil_data$camada_id == 2, "observacao_id"], layer.name = glue("20-40 cm<br>n = {n[2]}"))O próximo passo consiste em amostrar os valores das covariáveis nos locais de observação do solo. Para isso usamos raster::extract. Em seguida, especificamos as covariáveis de tipo categórico. Nesse estudo, temos apenas uma covariável categórica, especificamente, aquela que identifica os piquetes de pastoreio.
soil_data %<>%
cbind(., extract(covar_data, .)) %<>%
mutate(fields = as.factor(fields)) %T>%
print()## Warning in merge.data.frame(as.data.frame(x), as.data.frame(y), ...):
## column name 'x' is duplicated in the result## Simple feature collection with 627 features and 68 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 404071 ymin: 7606493 xmax: 406374 ymax: 7607536
## epsg (SRID): 32722
## proj4string: +proj=utm +zone=22 +south +datum=WGS84 +units=m +no_defs
## First 10 features:
## observacao_id observacao_data coord_sistema coord_precisao coord_fonte
## 1 8 xx-xx-2018 EPSG:32722 NA GPS
## 2 8 xx-xx-2018 EPSG:32722 NA GPS
## 3 35 xx-xx-2018 EPSG:32722 NA GPS
## 4 35 xx-xx-2018 EPSG:32722 NA GPS
## 5 48 xx-xx-2018 EPSG:32722 NA GPS
## 6 48 xx-xx-2018 EPSG:32722 NA GPS
## 7 54 xx-xx-2018 EPSG:32722 NA GPS
## 8 54 xx-xx-2018 EPSG:32722 NA GPS
## 9 57 xx-xx-2018 EPSG:32722 NA GPS
## 10 65 xx-xx-2018 EPSG:32722 NA GPS
## pais_id estado_id municipio_id amostra_tipo amostra_quanti amostra_area
## 1 BR SP Caiuá COMPOSTA 4 78.54
## 2 BR SP Caiuá COMPOSTA 4 78.54
## 3 BR SP Caiuá COMPOSTA 4 78.54
## 4 BR SP Caiuá COMPOSTA 4 78.54
## 5 BR SP Caiuá COMPOSTA 4 78.54
## 6 BR SP Caiuá COMPOSTA 4 78.54
## 7 BR SP Caiuá COMPOSTA 4 78.54
## 8 BR SP Caiuá COMPOSTA 4 78.54
## 9 BR SP Caiuá COMPOSTA 4 78.54
## 10 BR SP Caiuá COMPOSTA 4 78.54
## camada_id amostra_id camada_nome profund_sup profund_inf
## 1 1 277648/2018 - 0 20
## 2 2 277649/2018 - 20 40
## 3 1 277650/2018 - 0 20
## 4 2 277651/2018 - 20 40
## 5 1 277652/2018 - 0 20
## 6 2 277653/2018 - 20 40
## 7 1 277654/2018 - 0 20
## 8 2 277655/2018 - 20 40
## 9 1 277656/2018 - 0 20
## 10 1 277657/2018 - 0 20
## fosforo_resina matorg_cromo carbono_xxx_xxx_xxx ph_cacl2 ph_smp
## 1 4 14 8 4.9 6.95
## 2 2 10 6 4.6 6.95
## 3 7 16 9 5.1 7.04
## 4 3 11 6 5.0 6.99
## 5 5 15 9 5.0 6.90
## 6 2 12 7 5.1 7.17
## 7 4 16 9 5.0 7.03
## 8 2 11 6 4.9 7.11
## 9 4 14 8 5.0 6.98
## 10 6 12 7 5.2 7.21
## potassio_resina calcio_resina magnesio_resina sodio_resina
## 1 1.1 13 8 0.2
## 2 0.6 15 6 0.3
## 3 2.7 22 9 0.3
## 4 0.7 25 6 0.2
## 5 0.8 23 8 0.1
## 6 0.6 26 6 0.2
## 7 1.0 23 8 0.2
## 8 0.9 18 7 0.1
## 9 1.2 14 9 0.2
## 10 2.8 15 9 0.3
## acidez_xxx_xxx aluminio_xxx_xxx indefinido ctc_soma_calc
## 1 16 0 16 38.3
## 2 16 0 16 37.9
## 3 14 0 14 48.0
## 4 15 0 15 46.9
## 5 16 0 16 47.9
## 6 12 0 12 44.8
## 7 14 0 14 46.2
## 8 13 0 13 39.0
## 9 15 0 15 39.4
## 10 12 0 12 39.1
## bases_soma_calc bases_saturacao_calc aluminio_saturacao_calc
## 1 22.3 58 0
## 2 21.9 58 0
## 3 34.0 71 0
## 4 31.9 68 0
## 5 31.9 67 0
## 6 32.8 73 0
## 7 32.2 70 0
## 8 26.0 67 0
## 9 24.4 62 0
## 10 27.1 69 0
## enxofre_fosfca boro_h2o cobre_dtpa_eaa ferro_dtpa_eaa manganes_dtpa_eaa
## 1 12 0.50 1.1 7 16.9
## 2 18 0.85 1.5 11 12.4
## 3 13 0.79 1.7 12 31.9
## 4 16 0.56 2.4 10 30.3
## 5 19 0.65 1.1 23 23.0
## 6 17 0.68 1.4 10 18.0
## 7 12 0.79 1.1 18 23.1
## 8 16 0.67 1.1 11 24.3
## 9 15 0.59 1.6 17 24.4
## 10 17 0.64 1.6 9 20.1
## zinco_dtpa_eaa indefinido_1 indefinido_2 indefinido_3 indefinido_4
## 1 0.8 2.9 33.9 20.9 0.5
## 2 0.4 1.6 39.6 15.8 0.8
## 3 2.0 5.6 45.8 18.8 0.6
## 4 2.4 1.5 53.3 12.8 0.4
## 5 1.2 1.7 48.0 16.7 0.2
## 6 0.7 1.3 58.0 13.4 0.4
## 7 2.1 2.2 49.8 17.3 0.4
## 8 1.4 2.3 46.2 17.9 0.3
## 9 2.9 3.0 35.5 22.8 0.5
## 10 1.8 7.2 38.4 23.0 0.8
## indefinido_5 indefinido_6 indefinido_7 indefinido_8 indefinido_9
## 1 0 41.8 11.8 1.6 7.3
## 2 0 42.2 25.0 2.5 10.0
## 3 0 29.2 8.1 2.4 3.3
## 4 0 32.0 35.7 4.2 8.6
## 5 0 33.4 28.8 2.9 10.0
## 6 0 26.8 43.3 4.3 10.0
## 7 0 30.3 23.0 2.9 8.0
## 8 0 33.3 20.0 2.6 7.8
## 9 0 38.1 11.7 1.6 7.5
## 10 0 30.7 5.4 1.7 3.2
## argila_xxx_xxx_xxx silte_xxx_xxx_xxx areia_xxx_xxx_xxx slope
## 1 150 121 729 9.333595
## 2 199 94 707 9.333595
## 3 136 116 748 2.196754
## 4 164 77 759 2.196754
## 5 143 63 794 8.461255
## 6 179 96 725 8.461255
## 7 106 130 764 6.504688
## 8 148 85 767 6.504688
## 9 135 73 792 4.335512
## 10 111 76 813 3.765682
## general_curvature plan_curvature profile_curvature twi
## 1 1.191522e-03 0.0000000000 5.880601e-04 8.803287
## 2 1.191522e-03 0.0000000000 5.880601e-04 8.803287
## 3 8.115219e-04 -0.0118819950 6.662967e-04 11.429557
## 4 8.115219e-04 -0.0118819950 6.662967e-04 11.429557
## 5 -3.963677e-05 -0.0002342223 -2.306561e-10 8.512419
## 6 -3.963677e-05 -0.0002342223 -2.306561e-10 8.512419
## 7 5.171993e-03 0.0066209063 2.141720e-03 10.438464
## 8 5.171993e-03 0.0066209063 2.141720e-03 10.438464
## 9 7.391510e-04 0.0124076148 -1.678845e-04 11.107411
## 10 -4.687386e-04 0.0015604370 -2.925080e-04 11.383060
## catchment_area catchment_slope elevation CE37 CE75
## 1 7.860869 8.689552 364.060 -26.73634 -12.25075
## 2 7.860869 8.689552 364.060 -26.73634 -12.25075
## 3 8.963313 5.121825 353.627 -29.43822 -20.60455
## 4 8.963313 5.121825 353.627 -29.43822 -20.60455
## 5 6.759497 8.442846 366.538 -24.67430 -12.50217
## 6 6.759497 8.442846 366.538 -24.67430 -12.50217
## 7 8.377203 6.523236 359.246 -28.00345 -11.50673
## 8 8.377203 6.523236 359.246 -28.00345 -11.50673
## 9 8.560753 6.270308 355.107 -27.95664 -17.69285
## 10 8.809062 5.174417 348.728 -30.59239 -19.07843
## ndvi_2013_04_22 ndvi_2016_01_25 ndvi_2018_04_20 ndvi_variation fields
## 1 0.6467864 0.6196061 0.5970826 0.02488825 1
## 2 0.6467864 0.6196061 0.5970826 0.02488825 1
## 3 0.6139395 0.6013683 0.5501115 0.03381151 2
## 4 0.6139395 0.6013683 0.5501115 0.03381151 2
## 5 0.6353082 0.5886072 0.5792264 0.03003930 1
## 6 0.6353082 0.5886072 0.5792264 0.03003930 1
## 7 0.6582047 0.6088807 0.5923355 0.03426698 1
## 8 0.6582047 0.6088807 0.5923355 0.03426698 1
## 9 0.6111407 0.6205299 0.5571111 0.03422792 2
## 10 0.6352875 0.5729970 0.5438666 0.04670201 3
## water_tank_dist x y geometry
## 1 521.6222 404353.5 7607527 POINT (404358.2 7607536)
## 2 521.6222 404353.5 7607527 POINT (404358.2 7607536)
## 3 433.1423 404558.5 7607498 POINT (404568.9 7607508)
## 4 433.1423 404558.5 7607498 POINT (404568.9 7607508)
## 5 511.0452 404324.2 7607498 POINT (404310 7607486)
## 6 511.0452 404324.2 7607498 POINT (404310 7607486)
## 7 458.9842 404441.4 7607498 POINT (404433.9 7607499)
## 8 458.9842 404441.4 7607498 POINT (404433.9 7607499)
## 9 442.3896 404500.0 7607498 POINT (404495.3 7607501)
## 10 433.9739 404646.4 7607498 POINT (404654.7 7607490)1.2.1 Análise exploratória espacial
soil_var <- "potassio_resina"Vamos fazer uma rápida análise exploratória espacial dos dados de uma das propriedades do solo a serem mapeadas, por exemplo, a saturação por bases, utilizada no cálculo da necessidade de calagem. Para isso, usamos a função pedometrics::plotESDA, que gera um histograma da distribuição dos dados, um gráfico de bolhas, um semivariograma experimental, e um mapa do semivariograma experimental. Usaremos apenas os dados da camada de 0-20 cm de profundidade.
tmp <-
soil_data %>%
dplyr::filter(camada_id == 1) %>%
cbind(., st_coordinates(.))
plotESDA(z = tmp[[soil_var]], lat = tmp[["Y"]], lon = tmp[["X"]])
tmp <-
soil_data %>%
dplyr::filter(camada_id == 2) %>%
cbind(., st_coordinates(.))
plotESDA(z = tmp[[soil_var]], lat = tmp[["Y"]], lon = tmp[["X"]])
1.3 Calibração dos modelos
O primeiro passo da calibração de um modelo preditivo consiste na definição de uma fórmula que represete a relação entre a variável depedente e as variáveis preditoras. Nesse caso, a variável dependente é uma propriedade do solo, enquanto as variáveis preditoras são as covariáveis ambientais. Além dos atributos de terreno, das informações sobre o uso da terra e as práticas de manejo e conservação do solo, das imagens de satélite e índices de vegetação, e das informações indiretas sobre as propriedades do solo, usaremos como covariáveis ambas as coordenadas espaciais e o código de identificação da camada amostrada. A última é uma covariável importante pois, como sabemos, os dados provém das camadas de 0-10 e 20-40 cm de profundidade. Em geral, a camada superficial apresenta níveis superiores de nutrientes. Assim, os modelos poderão utilizar a a covariável camada_id para explicar essas diferenças e fazer predições específicas para uma ou outra camada.
Ao total, são 19 covariáveis, sendo que a covariável dos piquetes de pastoreio (fields) se desdobra em tantas covariáveis quantos foram os seus níveis, menos 1. Assim, na prática, são cerca de 50 covariáveis.
f <- glue("{soil_var} ~ x + y + slope + plan_curvature + profile_curvature + general_curvature + twi + catchment_area + catchment_slope + elevation + CE37 + CE75 + ndvi_2013_04_22 + ndvi_2016_01_25 + ndvi_2018_04_20 + ndvi_variation + fields + water_tank_dist + camada_id") %>% as.formula()
f_idw <- glue("{soil_var} ~ 1") %>% as.formula()1.3.1 Regressão linear múltipla
fit_glm <-
soil_data %>%
train(
f, data = ., method = "glm",
trControl = trainControl(method = "LOOCV"),
na.action = na.omit) %T>%
print()## Warning in predict.lm(object, newdata, se.fit, scale = 1, type =
## ifelse(type == : prediction from a rank-deficient fit may be misleading## Generalized Linear Model
##
## 627 samples
## 20 predictor
##
## No pre-processing
## Resampling: Leave-One-Out Cross-Validation
## Summary of sample sizes: 626, 626, 626, 626, 626, 626, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 0.8017863 0.1850615 0.5751148summary(fit_glm)##
## Call:
## NULL
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5808 -0.4298 -0.1492 0.2935 4.4359
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.696e+03 2.511e+03 -1.074 0.283297
## x 5.605e-04 7.876e-04 0.712 0.476963
## y 3.248e-04 3.369e-04 0.964 0.335411
## slope -8.704e-03 3.018e-02 -0.288 0.773110
## plan_curvature -1.668e+00 2.331e+00 -0.715 0.474615
## profile_curvature 2.601e+02 1.278e+02 2.036 0.042217 *
## general_curvature -6.682e+01 5.343e+01 -1.251 0.211510
## twi -1.786e-01 1.435e-01 -1.244 0.213858
## catchment_area 2.844e-02 6.061e-02 0.469 0.639048
## catchment_slope -7.990e-02 6.148e-02 -1.299 0.194295
## elevation 7.575e-03 2.829e-02 0.268 0.788964
## CE37 -1.491e-03 1.230e-02 -0.121 0.903577
## CE75 -8.952e-03 2.055e-02 -0.436 0.663329
## ndvi_2013_04_22 5.945e+00 2.373e+00 2.505 0.012515 *
## ndvi_2016_01_25 -2.299e+00 2.084e+00 -1.103 0.270469
## ndvi_2018_04_20 -3.175e+00 2.764e+00 -1.149 0.251190
## ndvi_variation -8.297e+00 5.681e+00 -1.461 0.144694
## fields2 -6.132e-03 3.164e-01 -0.019 0.984544
## fields3 1.690e-01 3.833e-01 0.441 0.659441
## fields4 -9.128e-02 5.296e-01 -0.172 0.863222
## fields5 -1.407e-01 6.091e-01 -0.231 0.817404
## fields6 -1.249e-01 6.756e-01 -0.185 0.853347
## fields7 5.819e-02 7.350e-01 0.079 0.936921
## fields8 2.089e-01 8.338e-01 0.251 0.802283
## fields9 -2.047e-01 9.356e-01 -0.219 0.826870
## fields10 -1.264e-01 1.164e+00 -0.109 0.913554
## fields11 6.181e-01 1.243e+00 0.497 0.619283
## fields12 -2.297e-01 1.298e+00 -0.177 0.859579
## fields13 -8.374e-01 1.264e+00 -0.663 0.507904
## fields14 2.620e-01 1.169e+00 0.224 0.822762
## fields15 -4.517e-01 1.306e+00 -0.346 0.729487
## fields16 -8.241e-01 1.242e+00 -0.663 0.507353
## fields17 -4.609e-01 1.149e+00 -0.401 0.688376
## fields18 2.288e-01 8.780e-01 0.261 0.794494
## fields19 2.807e-01 7.703e-01 0.364 0.715680
## fields20 1.075e-01 7.071e-01 0.152 0.879229
## fields21 1.534e-02 6.753e-01 0.023 0.981880
## fields22 -3.950e-01 5.846e-01 -0.676 0.499514
## fields23 -8.331e-02 4.325e-01 -0.193 0.847338
## fields24 9.473e-02 3.618e-01 0.262 0.793552
## fields25 -2.302e-01 2.975e-01 -0.774 0.439228
## fields26 3.464e-01 1.325e+00 0.261 0.793809
## fields27 1.370e+00 8.523e-01 1.608 0.108468
## fields28 2.039e+00 1.060e+00 1.923 0.054924 .
## fields29 8.863e-01 9.811e-01 0.903 0.366726
## fields30 6.760e-01 7.048e-01 0.959 0.337887
## water_tank_dist -1.054e-03 3.099e-04 -3.401 0.000718 ***
## camada_id -5.391e-01 6.282e-02 -8.581 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.5646149)
##
## Null deviance: 482.91 on 626 degrees of freedom
## Residual deviance: 326.91 on 579 degrees of freedom
## AIC: 1469
##
## Number of Fisher Scoring iterations: 21.3.2 Floresta aleatória
Os parâmetros do modelo de floresta aleatória serão mantidos constantes. São eles:
splitrule = "variance"min.node.size = 5mtry = floor(sqrt(p))
p <- fit_glm$finalModel %>% coef() %>% length() - 1
fit_rf <-
soil_data %>%
as.data.frame() %>%
train(
f, data = ., method = "ranger", importance = "impurity", num.trees = 500,
tuneGrid = data.frame(mtry = floor(sqrt(p)), splitrule = "variance", min.node.size = 5),
trControl = trainControl(method = "LOOCV"), na.action = na.omit) %T>%
print()## Random Forest
##
## 627 samples
## 19 predictor
##
## No pre-processing
## Resampling: Leave-One-Out Cross-Validation
## Summary of sample sizes: 626, 626, 626, 626, 626, 626, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 0.7397384 0.3040789 0.5375769
##
## Tuning parameter 'mtry' was held constant at a value of 6
## Tuning
## parameter 'splitrule' was held constant at a value of variance
##
## Tuning parameter 'min.node.size' was held constant at a value of 5fit_rf$finalModel## Ranger result
##
## Call:
## ranger::ranger(dependent.variable.name = ".outcome", data = x, mtry = min(param$mtry, ncol(x)), min.node.size = param$min.node.size, splitrule = as.character(param$splitrule), write.forest = TRUE, probability = classProbs, ...)
##
## Type: Regression
## Number of trees: 500
## Sample size: 627
## Number of independent variables: 47
## Mtry: 6
## Target node size: 5
## Variable importance mode: impurity
## Splitrule: variance
## OOB prediction error (MSE): 0.5493045
## R squared (OOB): 0.2879363par(mar = c(2, 8, 2, 2) + 0.1)
ranger::importance(fit_rf$finalModel) %>%
sort() %>%
barplot(horiz = TRUE, las = 1, col = "firebrick", cex.names = 0.75)
1.3.3 Interpolação ponderada pelo inverso da distância
# fit_idw <-
# soil_data %>%
# cbind(., st_coordinates(.)) %>%
# as.data.frame() %>%
# dplyr::filter(camada_id == 1) %>%
# gstat(id = "soil", formula = f_idw, locations = ~ X + Y, data = ., nmax = 7, set = list(idp = .5)) %T>%
# print()
# cv_idw <-
# soil_data %>%
# cbind(., st_coordinates(.)) %>%
# as.data.frame() %>%
# dplyr::filter(camada_id == 1) %>%
# gstat.cv(fit_idw, formula = f_idw, locations = ~ X + Y, data = ., nmax = 7, set = list(idp = .5),
# verbose = FALSE, debug.level = 0)
fit_idw1 <-
soil_data %>%
cbind(., st_coordinates(.)) %>%
as.data.frame() %>%
dplyr::filter(camada_id == 1) %>%
gstat(id = "soil", formula = f_idw, locations = ~ X + Y, data = ., nmax = 7, set = list(idp = .5)) %T>%
print()## data:
## soil : formula = potassio_resina`~`1 ; data dim = 399 x 69 nmax = 7
## set idp = 0.5;
## ~X + Y
## <environment: 0xd7e1548>cv_idw1 <-
soil_data %>%
cbind(., st_coordinates(.)) %>%
as.data.frame() %>%
dplyr::filter(camada_id == 1) %>%
gstat.cv(fit_idw1, formula = f_idw, locations = ~ X + Y, data = ., nmax = 7, set = list(idp = .5),
verbose = FALSE, debug.level = 0)
# sapply(seq(2, 20, by = 4), function (nmax) {
# soil_data %>%
# cbind(., st_coordinates(.)) %>%
# as.data.frame() %>%
# dplyr::filter(camada_id == 1) %>%
# krige.cv(formula = f_idw, locations = ~ X + Y, data = ., nmax = nmax, set = list(idp = .5),
# debug.level = 0) %>%
# dplyr::select(residual) %>%
# mutate(residual = residual^2) %>%
# summarize(rmse = sqrt(mean(residual))) %>%
# unlist()
# })
fit_idw2 <-
soil_data %>%
cbind(., st_coordinates(.)) %>%
as.data.frame() %>%
dplyr::filter(camada_id == 2) %>%
gstat(id = "soil", formula = f_idw, locations = ~ X + Y, data = ., nmax = 7, set = list(idp = .5)) %T>%
print()## data:
## soil : formula = potassio_resina`~`1 ; data dim = 228 x 69 nmax = 7
## set idp = 0.5;
## ~X + Y
## <environment: 0xe910400>cv_idw2 <-
soil_data %>%
cbind(., st_coordinates(.)) %>%
as.data.frame() %>%
dplyr::filter(camada_id == 2) %>%
gstat.cv(fit_idw2, formula = f_idw, locations = ~ X + Y, data = ., nmax = 7, set = list(idp = .5),
verbose = FALSE, debug.level = 0)1.4 Validação dos modelos
Estatísticas do erro de predição na validação cruzada, levando em consideração apenas as observções da camada de 0-20 cm de profundidade. As estatísticas do erro são erro médio (ME), erro mediano (MedE), erro quadrático médio (MSE), erro absoluto médio (MAE), quantidade de variância explicada (AVE), e quantidade de variância explicada modificada (mAVE). A quantidade de variância explicada (amount of variance explained) é equivalente ao coefficiente de deteminação. Já mAVE consiste no uso de erros e resíduos absolutos, ao invés do seu quadrado, para calcular a quantidade de variância explicada. Trata-se de uma medida robusta à presença de observações discrepantes que inflam o AVE. Ambas as medidas representam o quão eficiente o modelo é em explicar a variância dos dados em comparação ao uso da média dos dados.
soil_data %>%
mutate(idw = c(cv_idw1$soil.pred, cv_idw2$soil.pred) - c(cv_idw1$observed, cv_idw2$observed),
glm = fit_glm$pred$pred - .[[soil_var]],
rf = fit_rf$pred$pred - .[[soil_var]]) %>%
as.data.frame() %>%
dplyr::select(idw, glm, rf) %>%
stack() %>%
mutate(residual = rep(mean(soil_data[[soil_var]]) - soil_data[[soil_var]], 3)) %>%
dplyr::group_by(ind) %>%
dplyr::summarise(
ME = mean(values),
MedE = median(values),
MSE = mean(values^2),
MAE = mean(abs(values)),
AVE = 1 - sum(values^2) / sum(residual^2),
mAVE = 1 - sum(abs(values)) / sum(abs(residual)))## # A tibble: 3 x 7
## ind ME MedE MSE MAE AVE mAVE
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 idw -0.000240 0.160 0.621 0.570 0.194 0.124
## 2 glm -0.00460 0.160 0.643 0.575 0.165 0.116
## 3 rf 0.00987 0.223 0.547 0.538 0.290 0.1741.5 Predição espacial
covar_data$camada_id <- 1pred_glm <- covar_data %>% predict(fit_glm, factors = list(fields = levels(soil_data$fields)))## Warning in .local(object, ...): not sure if the correct factor levels are
## used herepred_rf <- covar_data %>% predict(fit_rf, factors = list(fields = levels(soil_data$fields)))## Warning in .local(object, ...): not sure if the correct factor levels are
## used herepred_idw <- interpolate(covar_data, fit_idw1, xyNames = c("X", "Y"))## [inverse distance weighted interpolation]par(mfrow = c(2, 2))
plot(pred_glm, main = glue("{soil_var} - glm"))
plot(pred_rf, main = glue("{soil_var} - rf"))
plot(pred_idw, main = glue("{soil_var} - iwd"))
mapview(pred_glm, alpha.regions = 1) +
mapview(pred_rf, alpha.regions = 1) +
mapview(pred_idw, alpha.regions = 1) +
soil_data %>% dplyr::filter(camada_id == 1) %>% dplyr::select(soil_var) %>% mapview() +
mapview(contour_lines)tmp <-
soil_data %>%
cbind(., st_coordinates(.)) %>%
mutate(
res_rf = fit_rf$finalModel$predictions - .[[soil_var]],
res_glm = fit_glm$finalModel$fitted.values - .[[soil_var]]) %>%
dplyr::filter(camada_id == 1)
plotESDA(z = tmp[["res_glm"]], lat = tmp[["Y"]], lon = tmp[["X"]])
plotESDA(z = tmp[["res_rf"]], lat = tmp[["Y"]], lon = tmp[["X"]])
tmp <-
soil_data %>%
cbind(., st_coordinates(.)) %>%
mutate(
res_rf = fit_rf$finalModel$predictions - .[[soil_var]],
res_glm = fit_glm$finalModel$fitted.values - .[[soil_var]]) %>%
dplyr::filter(camada_id == 2)
plotESDA(z = tmp[["res_glm"]], lat = tmp[["Y"]], lon = tmp[["X"]])
plotESDA(z = tmp[["res_rf"]], lat = tmp[["Y"]], lon = tmp[["X"]])